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Abstract 

We analyze stability of a system which contains an harmonic os- 
cillator non-linearly coupled to its second harmonic, in the presence 
of a driving force. It is found that there always exists a critical ampli- 
tude of the driving force above which a loss of stability appears. The 
dependence of the critical input power on the physical parameters is 
analyzed. For a driving force with higher amplitude chaotic behav- 
ior is observed. Generalization to interactions which include higher 
modes is discussed. 
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1 Introduction 



In a series of experiments the motion of the surface of a superfluid liquid in 
a cylindrical vessel was studied. This motion was induced by standing waves 
of second sound propagating in the bulk of the liquid. Above a critical value 
of the input power the motion has lost stability ]5j, §. 

To account for this loss of stability we analyzed a model that explained 
this phenomenon ||, and found it to be in a good agreement with the ex- 
perimental results. The model is general enough to account for the loss of 
stability in other wave systems. 

2 The Model 

The model consists of two non-linearly coupled harmonic oscillators, of which 
one is coupled to an external driving force. First we would like to justify the 
use of two oscillators, with frequencies close to uj and 2o>, for describing the 
physics of systems such as the one above (u would be the frequency of the 
driving force). We assume that in the linear approximation the free, non- 
dissipative (classical) theory is given by the Hamiltonian: 

oo 

H = u n a* n a n (1) 

n=l 

where a n is the (complex) amplitude of the n th mode, and a* is its complex 
conjugate. Dissipation and driving force would be added in the following. 
The modes are the eigenfunctions of the wave equation with the appropriate 
Sturm-Liouville boundary conditions. 

We neglect terms higher then cubic from the Hamiltonian, as well as terms 
which are far from resonance, and therefore have small coupling constant 
The Hamiltonian turns to: 

oo oo 

H = Y1 w na* n a n + (^k,i; m akaia* m + c.c.) (2) 

n=l k + l — m~0 

k,l,m=l 

where c.c. stands for complex conjugate, and Xk,i ]m = \,k;m- We would like 
now to couple an external driving force to one of the modes. We keep in 
mind that, in order to describe a physical problem, attenuation should be 
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added as well. The modes which are not strongly coupled to the excited mode 
would decay. Again we assume that, for describing the onset of instability 
a minimal number of modes is needed. Therefore we take the excited mode 
and the mode with frequency which is closest to twice the frequency of the 
first one. With the harmonic driving force the Hamiltonian takes the form: 

H = u d a* d a d + uj 2d a* 2d a 2d + (Xa d a* 2d + c.c.) + (fe luJt a* d + c.c.) (3) 

where u is the frequency of the driving force, which should be close to u> d in 
order to resonate it. 
We use: 

ia d = — (4) 
da* d 

which is Hamilton's equations in the amplitude formalism |l|], to derive 
the equations of motion: 



id d = u d a d + 2X*a* d a 2d + fe tujt (5) 
ia 2d = u 2d a 2d + Xa 2 d (6) 

The equations are invarant under the transformation: 

a d -> a d e^ +e ) 
a 2d -> a 2d e^-V 

f f e i(<P+0) 

It is therefore possible to eliminate two independent phases from the 
equations, so we can choose A and / to be real. 

We add now dissipative terms to the equation in the usual maner . The 
equations now become: 



ia d = (u d -i ld )a d + 2\a* d a 2d + fe iu>t (8) 
ia 2d = (u 2d - i~f 2d )a 2d + Xa 2 d (9) 

where 7 are the dissipation constants. 
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The final stage before analyzing the equations is to introduce the "slow 
variables" to eliminate the time dependence. Under the transformation: 



a d -> a d e iuJt 

&2d — ► d 2d e 



the equations get the form: 



ia d = (A d - i-y d )a d + 2Xa* d a 2d + f (10) 
ia 2d = (A 2d - ij 2d )a 2d + Aa 2 , (11) 

where A d = uj d — uj and A 2d = u 2d — 2uj are the frequency detunings with 
respect to the driving force. 

We note that in deriving the Hamiltonian (|3]) we have neglected one non- 
linear term which is of the same order with respect to a d as the one that we 
have kept, namely na 2 d a* d . We would like to justify it. For systems in which 
K << — this term is small, but it turns out that even for k much larger 
the importance of this term is not crucial. Note, that uj is absent in the 
equations. From dimensional considerations k can appear at the equations 
only as kA, kj, this is the term that has to be of the order of A 2 . Hereafter we 
analyze three aspects of the model: stationary solutions, stability, numerical 
calculations. For the stationary solutions it is easy to verify that the effect of 
k is merely to renormalize A 2d , j 2d . This is the well known effect of shifting 
the resonance 10. We have seen that k is not of a big importance, even for 
k > , — , in the stability analysis as well as in our numerical calculations. 
We will not include this term in what follows. 

Although the model we use is a very simplified one, it still contains five 
parameters in addition to the driving force amplitude /. Not all the pa- 
rameters are important. The amplitude / of the driving force is an effective 
expression which in fact is a function of A d , moreover, the driving force cou- 
ples to all other modes as well, and we may neglect all other couplings only 
when the one that we are left with is the dominant one. For this to be the 
case we must have u ~ u d , that is, A d is small compared to all other pa- 
rameters with dimensions of frequency. The value of A 2d will be dictated 
by geometry. Both, our analytical, as well as numerical results depend on 
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this assumption. In most physical systems there is a relation between 7^ 
and 72d. We shall assume that this two parameters are of the same order of 
magnitude. 



3 Stationary solutions 

We begin our analysis by finding the fixed points of the equations, i.e., solving 
the equations: 

(A, - i ld )a d + 2Xa d a 2d + f = (12) 

(A 2d - i^2d)a 2 d + Aa 2 , = (13) 

We eliminate a 2d from the second equation, and substitute in the first one to 
get: 

(A rf - i^ d )(A 2d - ij 2d )a d - 2A 2 a (i |a ( i| 2 = -(A 2d - ij 2d )f (14) 
The equation for C = -, 2A ! A — rlaJ 2 turns now to: 



(((±l) 2 + (3)( = h (15) 



fes 2( 7 2 d + AL)A 2 
\7dl2d- A d A 2d | 3j 



where: 



is the scaled force, and 



The sign in equation flip]) coincides with the sign of 7^72^ — ArfA 2 d- 
This equation has either one or three solutions. For a given value of h 
the equation would have three solutions if and only if: 

ldl2d - A d A 2d < (18) 
0<P<\ (19) 
A[l + 9/3 - (1 - 3/3)1] < h < ^[1 + 9/3 + (1 - 3/3)i] (20) 
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In the following, it would be illustrated that, when three solutions are 
present, the middle one is non-stable, as may be expected. 

We note that the situation of three solutions is, in a sense, non-physical. 
Id, l2d are positive, we can therefore use (pT7|)- (pT9|) to deduce: 




But this suggests that ^ d < A^, which contradicts our assumptions. In 
this region of parameters our model is not adequate. 



4 Stability 

To check whether the stationary solutions are stable we linearize the equa- 
tions around these solutions, and check whether small perturbations grow 
or decay. To simplify the calculations, we recall the symmetry (|7]) and use 
it with + 39 = to redefine the stationary value of the first mode, a d \ 
to be real, without altering A. The change of / is not important since / 
will be absent from the linearized equations. We substitute in the linearized 
equations: 



■2 



\ (o) 

"2 = -rV < 22) 

The stability is now determined by a d °\ Also, to simplify the notations, 
we will use a d , a 2 d rather then Sa d , 5a 2 d for the deviations from the stationary 
solution. 

The linearized equations are: 



ia d = {A d - i ld )a d + 2A(4 0) a 2d - - Xa ' d , a*) (23) 

A 2 d - ll2d 

ia 2 d = (A 2 d - ij2d)a2d + 2Xa d 0) a d (24) 

multiplying by —i and separating to real and imaginary parts, we get the 
differential equation: 



(O) 2 
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/ Re(a d ) \ 
d Im{a d ) 
(It Re(a 2d ) 

\Im(a 2d ) J 

( ~ld-Pl2d 

-A d + pA 2d 




V 



-2afX 



2af\* 



A d + pA 2d 

~ld + Vl2d 



4 0) 





2af\ 



2ai 0) A\ / Re (a d ) \ 



~l2d 

-A 2d 





A 2d 

~l2d ) 



Im(a d 
Re(a 2d ) 
\lm{a 2d )) 



(25) 



where p — ■>,,■.. 

To ensure stability we shall require that the real part of all the eigen- 
values of this matrix is negative. We find the coefficients of the characteristic 
polynomial u 4 + au 3 + bu 2 + cu + d, to be: 



2(7rf + l2d) 
4A 4 

'lld + ^ld 

8A 4 72rf 



,(o) z 



+ 8A 2 ai 0) ' 



(26) 

+ (id + Ad + lid + A ld + ^ldl2d) (27) 



„,2 _l A2~ a d -T as\ y] d -r j2dju d 
l 2d ~r >-± 2d 

+2[(ld + & 2 d )l2d + (l 2d + A 2 2d ) ld 



+ 8\ z {ld + l2 d )a ( d y + 



(28) 



d = 12A 4 4 0) - + 8A 2 ( 7(i 72d-A d A 2d )ar" + (7^ + A 2 )( 7 ^ + Ay (29) 

To ensure that all the roots of this polynomial have negative real part we 
use the Hurwith-Routh criterion 01(11: 



(0)' 



a > 
b > 
d > 









abc > c + a d 



(30) 
(31) 
(32) 
(33) 



The condition ([30]) is trivial for a physical problem. The condition (|31|) 

(o) ^ 

is a quadratic equation with respect to a d , and is easily solved to give: 
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a (0)2 < 7| 1 +AL (4 + 



4*-^\ 5+ iid + a 2 } (3) 



i 2d 

The third condition, fl32|) , is again a quadratic equation with respect to 

, but with positive, rather then negative coefficient of a\. 
It is easily seen that for a negative d to occur at the physical region: 

(q\2 

a d > 0, we need to have: 

7d72d > A d A 2d (35) 
When this condition is fulfilled, an unstable region appears when: 

P < I (36) 

Direct solution of the quadratic equation then shows that the central 
region of solutions coincides exactly with this unstable region (p0|). As men- 
tioned above, this region is not physically important. 

We combine ©-(gD and ([§), and define: 

z = A 2 af 2 (37) 

to get the last inequality: 

a z A + aiz 3 + a 2 z 2 + a 3 z + a 4 > (38) 

where: 



ao " did + A L) 2 (39) 



64( 7d + l2d f 

ai = V +A 2 (40) 

12d + ^2d 



a 2 = 2 l a " a 2 (A 2 d - A 2 - ( 7d + 7 2,) 2 ) (41) 

72d ' Id 

a 3 = 16( 7d + 7 2d) 2 [(7d + 72d) 2 + (A d + A 2d ) 2 )] (42) 
a 4 = ^ldl2d[(ld + l2d) 2 + (A d + A M ) 2 ][( 7d + 7m ) 2 + (A d - A M ) 2 )(]43) 
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It is seen that for all parameter values there exists an open neighborhood 
of zero in which the stationary solution is stable. It is very tedious to solve 
the inequality for the general case. We solve it for two special cases, the 
one-dimensional geometry, and the cylindrical wave with a large Q-factor. 
Both with reflecting boundary conditions. 

We remind the assumption: A d << j d . It is natural to assume that 7^ 
and j 2 d are of the same order of magnitude. In a wide class of cases 7 oc u 2 , 
and therefore: 

J2d ^ 47 d (44) 

We shall consider this case for both geometries. The value of A2d is dictated 
by geometry. 

For the one dimensional geometry the d th mode is cos(d^), where L is 
the length of the vessel. This dependence gives: 

A 2 d = u 2d -2u = 2A d - (2u d - u 2d ) = 2A d - c(2k d - k 2d ) = 

= 2A d - c(2f - f*) = 2A d l4Dj 

We therefore, have for the one-dimensional case: 

A d , A 2d « 7d ,7 2d (46) 

We define: 



x = 4 (47) 

Id 

l2d 



Id 

We use (PjP to get the inequality: 



(48) 



x 4 - 8(1 + s) 2 x 3 - \s 2 (l + sfx 2 + \s\l + sfx + -U 4 (l + sf > (49) 
2 4 16 



The solution of this inequality combined with (34) gives the final result: 

x< l -{s + s 2 ) (50) 
from which one easily finds an expression for the critical input power: 
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A 

or using (fl4|): 



/c = x \/72d o ( 51 ) 



/c - f 7, 2 (52) 
we now substitute j d = au; 2 to get: 

/ 22a2 ^ ^ 

fe ^ ( 53 ) 

A full description of the loss of stability for the specific problem can be 
obtained if we take into account the dependence of a and A on the relevant 
physical parameters, i.e. temperature. 

For a cylindrical vessel of radius R the modes are given by J n (kr) cos(n6>), 
where J n is the n th Bessel function, and k — - where c is the wave velocity. 
The boundary conditions force the relation k n ^ m R = Xn,m where Xn,m is the 
m th zero of J' n {x)- For simplicity we shall consider here only the Jo modes. 

The value of A 2d is dictated by the Bessel asymptotics: 

Xm = Xo,m ^ nn + j (54) 

using which we get: 

A 2d = uj 2 d -2u = 2A d - (2u> d - u 2d ) = 2A d - c(2k d - k 2d ) ~ 



2A d - |(2xd - X2d) = 2A d - |(2(<k + f ) - (2dn + f )) = (55) 
2A rf - ^ ~ 2A 



4R — M+l 

Since cu^ = 2Q / y d , the higher is Q, the higher are the values of d for which 
the inequality 

A 2d » 7d (56) 

holds. 

We solve now equation ( |38|) for the case: 

A d « 7^,72^ << A 2d (57) 
We define s as before, but now: 

* = 7TT (58) 
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and get the equation: 



4 (l + s) 2 , 1 9 (1 + s) 2 1 
s 2 4s 16 

When this condition is combined with ([53]) we get: 

x<^{v + V2uv) (60) 

where w = , and v = u — \/u 2 — 4. We use QS4]) to get: a; < 0.59. For 
other values of s there are only small changes in the result. In all cases the 
critical value is in the range: 0.5 < Xq < 0.65. The maximal value is attained 
at s = 1, and the minima are at x = 0, x — > oo (note that Xq(s) = x (^)). 
The critical input power f c may be calculated now: 

/• = ^ = < 61 > 

A full description of the loss of stability in this geometry can be obtained 
if we take into account the dependence of A and c on the relevant physical 
parameters. 



5 Beyond - Numerical calculations. 

Some questions arise. Does the system always reach the stationary solution 
in the stable region? What happens above the stable region? In what way 
would the theory be modified if we include the full Hamiltonian ([D? 

We solved the equation numerically with parameters suitable to describe 
the cylindrical geometry: 

A d = A 2d = 1500 

Id = 30 7 2d = 120 (62) 
A = 5400 

with initial conditions: 

a d (t = 0) = a 2d (t = 0) = (63) 
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The results for other values of parameters may be very similar due to the 
scaling properties of the equation discussed above. 

For small enough values of / the system reaches the stationary solution 
after some travelling in phase space (Fig. |l|). For / ~ 0.3/ c with the initial 
conditions above, the system escapes the basin of attraction of the fixed 
point, and rather approaches a limit cycle (Fig. |2|). The basin of attraction 
of the stationary solution shrinks to zero as the instability is approached. 
This limit cycle is not unique. By choosing various initial conditions other 
limit cycles can be approached. In the higher / regime the behavior is harder 
to determine. 




Figure 1: The system approaches the fixed point for / = 50, the position of 
the fixed point is indicated. 

It is easy to prove that the motion of the system is bounded in its phase 
space, and that the volume in phase space decays exponentially with decay 
factor 2(7 d + 7 M ). 

A necessary condition for chaos to evolve is that the system will be unsta- 
ble locally. Our analysis shows that the phase of ad, a2d is irrelevant to this 
question. Given the value of the parameters, the potentially chaotic regions 
are defined in the (|arf| 2 , |«2d| 2 ) plane. Our calculations show that the region 

M 2 A 2 » 7 2 ,A 2 (64) 
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Figure 2: The system approaches a limit cycle for / = 100 which is below 
the critical value. 

is always locally unstable. The numerical calculations show that when / is 
increased the system enters this region, bifurcations appear, as in the usual 
root towards chaos. For large enough / chaos will evolve. 

In (Fig. |3|) we see the bifurcations for / = 500. 

Chaos evolves for / ~ 506 as we see in (Fig. ^). 

When more modes are added to the system the behavior changes. The 
projection of, say, the 3-mode system on the (4dim) phase space of 2 modes 
gives, in general, trajectories which are very different from the original ones. 
Yet, we argue that the main conclusion does not change. Indeed, if we exam- 
ine the original set of equations (]10|) , (|TTD we note that the transformation: 

a d -> aa d 

a 2d oia 2 d / CK \ 

/ - af < 65 > 
A ^ -A 

a 

which is a generalization of (|7|), leaves the equations invariant. We could 
deduce from here that f c oc jk From dimensional considerations / should 
be proportional to 7 2 , A 2 . It is seen that for the one dimensional case the 
leading behavior is: f c oc 7 2 , while for the large Q-factor case f c oc A% d . All 
our calculations were in fact needed just to illustrate that there is only one 
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transition from stability to instability (i.e. no instable windows), to validate 
the assumption that the largest constant with frequency dimensions is not 
absent from the expression for f c , and to calculate x . When we add new 
modes, new constants are added to the system. Since from ([54] ) we get that 
for all j Aj oc -|, these constants do not cause a problem. The same is true 
for the one-dimensional case. As for the new X's, if they scale in some way, 
e.g. if 

/ 777 

X k>l . m = f(T,R,...)h(-,j)k u (66) 

where f(T,R,...) is any function of all physical parameters, but the wave 
length, are constants, and u is an exponent, then the symmetry still 

holds and then, given that the general picture remains the same, all that we 
need to change is the value of x . This necessary modification of x , plus the 
shrinking of basin of attraction, which effectively lowers Xq, suggests that this 
part of our calculations is not reliable. Yet, the dependence of the critical 
input power on all physical parameters remains the same even for the full 
Hamiltonian (Q). There are values of u, h(j-, for which other predictions, 
such as the distribution of the chaotic regions of the 2-mode system would 
not be dramatically changed as well. More extensive investigation of this 
system is, therefore, highly desirable. 
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